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Introduction 


The  problem  of  the  formation,  growth,  and  dissipation  of  a  cumulus  cloud 
is  so  complex,  that  it  has  so  far  defied  solution.  The  cumulus  cloud  is  a 
three-dimensional  space  phenomenon,  and  should  be  treated  as  such.  'During 
the  last  few  years,  several  three-dimensional  convection  models  have  been 
constructed  CWilliams,  1969:  Deardorff,  1970:  Fox,  1972:  Steiner,  1973: 
Wilhelmson,  1974:  Miller  and  Pearce,  1974:  Schlesinger,  1975:).  The  number 
of  models  has  not  been  great,  mainly  because  they  demand  cumbersomely  large 
amounts  of  computing  storage  and  expenditures  if  adequate  spatial  resolution 
is  to  be  achieved.  Indeed,  even  in  the  most  sophisticated  model,  i.e.,  that 
of  Schlesinger  (1975),  the  calculations  had  to  be  limited  to  preliminary  “mid¬ 
get"  experiments  with  sufficiently  few  grid-points  to  circumvent  these  demands. 

By  contrast,  in  recent  years,  a  number  of  two-dimensional  models  have 
been  developed  for  the  study  of  deep  precipitating  convective  clouds  (Takeda, 
1966:  Orville  and  Sloan,  1970b:  Schlesinger,  1973a,  b:  Hane,  1973).  Two- 
dimensional  axisynmetric  models  have  also  been  developed  (Ogura,  1963:  Murray, 
1970:  Soong  and  Ogura,  1973).  These  models  demand  far  less  computer  time 
than  three-dimensional  models,  and  are  able  to  reproduce  many  important  features 
of  the  cumulus  cloud.  Even  simpler,  less  time  consuming  models  are  "quasi- 
two-dimensional"  models  (Weinstein  and  Davis,  1968:  Berkofsky,  1974:  Simpson 
and  Dennis,  1975).  Such  models  have  given  remarkably  reasonable  results  in 
cumulus  dynamics  studies,  especially  when  used  for  operational  purposes. 

In  view  of  the  above  comments,  and  in  view  of  the  fact  that  little  work 
has  beer,  done  in  this  area  relative  to  desert  cumulus,  the  approach  to  be  taken 
here  will  be  to  develop,  at  first,  a  "quasi -two-dimensional"  model,  and  then 
go  on  to  more  sophisticated  models. 

The  purpose  of  the  current  proposed  research  is  to  develop  a  model  for  the 
prediction  of  the  formation,  growth  and  dissipation  of  a  cumulus  cloud  in  a 
desert  region.  The  investigation  of  cumulus  clouds  is  important  because  such 
clouds  release  large  amounts  of  energy,  especially  in  tropical  and  sub-tropical 
latitudes,  and  thereby  profoundly  influence  the  atmospheric  circulation  in 
those  regions.  Furthermore,  there  are  certain  practical  advantages  in  being 
able  to  predict  the  development  of  desert  cumulus  clouds.  As  rainfall  is 
sparse  and  spotty  in  desert  areas,  there  would  accrue  profound  advantages  to 
the  agriculture  and  other  aspects  of  the  economy  of  such  regions  if  the 
occurrence  or  non-occurrence,  and  the  location  of  precipitation,  could  be  pre¬ 
dicted.  Additionally,  the  ability  to  predict  intensity  and  amount  of  rainfall 
would  be  extremely  important  in  relation  to  desert  floods.  Finally,  a  by¬ 
product  of  such  studies  would  be  the  ability  to  consider  weather  modification 
possibilities. 

Cloud  Model  and  Equations 

Ue  shall  consider  a  cloud  model  in  which  the  variables  are  functions  of 
height,  z,  and  time,  t.  The  model  is  axisymmetric,  but  the  radius  of  the  cloud 
is  itself  allowed  to  vary  with  height.  Thus  the  cloud  model  has  the  character¬ 
istics  of  a  "plume"  model.  The  system  to  be  studied  has  the  following  main 
assumptions  implicit  in  it. 
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1.  Cloud  is  condensed  water  that  fully  shares  the  air  motion. 

2.  Clou^orms  in  rising  saturated  air  and  evaporates  in  descending  air  at  the 

rate-wP^where  is  the  saturatinn  mixing  ratio  of  water  in  air,  is 

the  vertical  velocity  of  the  air,  f  is  the  air  density.  Cloud-containing  air 
is  always  saturated,  and  unsaturated  air  never  contains  clouds. 

3.  Cloud  changes  to  raindrops  that  are  distributed  in  size  according  to  an 
Inverse  exponential  distribution  at  the  rate-ii.fpOj-O  where  the  magnitude 

of  ki  and  a  may  be  selected  to  simulate  various  processes  and  rates.  Q'  is 
the  “cloud  water"  -  water  which  is  condensed  in  the  updraft,  forms  very'small 
droplets  which  have  a  negligible  terminal  velocity  and  are  thus  totally  carried 
along  by  the  updraft.  is  called  "conversion". 

4.  Precipitation  particles  once  formed  are  assumed  to  be  distributed  in  size 
according  to  an  inverse  exponential  law  and  to  collect  cloud  particles  or 
evaporate  in  subsaturated  air  according  to  approximations  to  the  natural 
accretion  and  evaporation  processes.  This  is  called  the  "collection"  process. 

To  derive  the  appropriate  equations  for  the  model,  we  begin  by  defining 

dmc  *  dm^  -  dm_  (1) 

[see  Ooyama  (1971)],  where  dm^.  represents  the  change  of  mass  of  a  rising  plume 
as  the  net  result  of  entrainment  dm+  and  detrainment  dm_.  Thus,  as  the  mass 
rises  from  z  to  z+dz  its  mas  increases  from  m  to  m-Klm.  In  general,  it  is 
assumed  that  entrained  air  brings  q/g^into  the  plume  (‘Vis  a  scalar,  and  e 
means  environment)  and  the  detrained  air  takes <v^;(c  means  cloud)  out  of  it. 
Therefore,  if  cyis  a  conservative  property,  the  general  conservation  law  is 
wri tten : 

<L(af^vv\c')  q«  <<-vvfc+  — 
or 

y^c  clQ'c  -<he) 

Writing  T  (temperature)  for  and  multiplying  Equation  (2)  by  Cp,  we  find 

-  Cp  VV\f  clT  =  p  ^  ^C4-  (3) 

where  {dh)^  means  sensible  heat  gain  or  loss  due  to  entrainment,  Cp  is  specific 
heat  at  constant  pressure.  Similarly,  ^ 

UvAc  ^  ^4) 

where  (dh)^  is  a  heat  loss,  L  is  latent  heat  of  condensation,  q  is  mixing  ratio. 
The  heat  loss  is  due  to  evaporation  of  liquid  water  necessary  to  resaturate 
the  air  which  has  become  slightly  subsaturated  due  to  entrainment. 

When  the  supercooled  water  in  the  parcel  is  frozen,  there  are  two  more  non 
adiabatic  heat  transfer  terms  to  be  considered.  The  first  is  the  heat  gained 
due  to  the  release  of  latent  heat  of  fusion. 


(5) 


where  Lf  is  the  latent  heat  of  fusion,  and  Q  is  the  amount  of  supercooled 
water  tnat  is  frozen,  in  units  of  grains  of  water  per  gram  of  air.  Here 

-  Qc 

where  Qc^  is  cloud  water,  i.e.  the  water  condensed  in  the  updraft,  having  , 
negligible  terminal  velocity,  and  totally  carried  along  by  the  updraft. 
is  called  hydrometeor  water,  i.e.  a  class  of  drops  which  have  grown  from 
Qr  due  to  contained  condensation,  which  have  important  terminal  velocity. 
The  process  by  which  hydrcmeteor  water  is  formed  is  called  conversion. 

Before  freezing  occurs,  the  liquid  water  and  water  vapor  in  the  parcel 
are  in  vapor  equilibrium  at  the  saturation  vapor  pressure  over  a  water 
surface.  Immediately  after  freezing  has  occurred,  a  similar  equilibrium 
must  be  achieved  between  the  newly  formed  ice  and  the  vapor.  Since  the 
saturation  vapor  pressure  over  an  ice  surface  is  less  than  over  a  water 
surface  at  the  same  temperature  and  pressure,  there  must  be  a  deposition 
of  vapor  onto  the  ice.  The  heat  gained  by  the  parcel  by  this  deposition  is 

(at),  ~  ''S  (  ice  (7) 

where  Ls  is  the  latent  heat  of  sublimation,  and  watervice  is  the 

difference  between  the  saturation  mixing  ratio  over  a  water  surface  and 
that  over  an  ice  surface. 

The  first  law  of  thermodynamics  is 

(AT - 

where  dh  is  the  heat  gained  or  lost  by  non-adiabatic  processes,  is  air 
density,  is  pressure. 

Using  the  hydrostatic  relation  only  at  this  point. 


(9) 

we  ran  substitute  Equations  (3),  (4),  (5),  (7),  (9)  into  Equation  (8)  to 
obtain 


^  - 


-UTe 


(10) 


V  ) 

where  Ai  is  the  heat  equivalent  of  work,  0  =  t  *  mechanical  equivalent  of 
heat,  e  =  0.621. 
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We  shall  derive  an  expression  for  u  later. 

The  above  Equation  (10)  for  the  rate  of  change  of  temperature  within  the 
cloud  reduces  to  the  equation  for  the  pseudo-adiabatic  lapse  rate  when  there 
is  no  entrainment  (v=0),  and  we  write 


_  _J_  dT 
^  d-i. 

where  w  is  the  vertical  velocity  of  the  parcel. 

The  well-known  Clapeyron  equation  may  be  written 


(12) 


+ 


€  t 
Tc 


(13) 


One  way  of  deriving  the  cloud  water  equation  is  to  apply  Equation  (1)  to  the 
sum  of  vapor  and  cloud  water  within  the  cloud  (where  <3,.  =  by  assumption). 


I 

Since  G<ce-®~  convective  cloud  water  in  the  environment. 


d  t 


iii 


(14) 


(15) 


If  we  assume  that  is  is  permissible  to  add  the  collection  and  conversion  terms 
(to  be  specified  explicitly  below),  even  though  thsy  destroy  the  conservation 
properties  to  some  extent,  we  get 


li 


--- 

-conversion  -  collection 


(16) 


A  more  sophisticated  treatment  starts  from  the  equations  of  messier  (1969). 


c  -  u  -V  ~  Ur  4-ur  (o,  \hl£.  - 

>7  ^  ^  '^T:  Si  ' 

4  Qc  UJ-  . C Qlc 


■Tit 


(17) 


and 


+  4,((5,-a)  C};^  (1 

V 

In  these  equations,  Qq,  Q|j,  Qg  are  the  cloud  water,  hydrometeor  water, 
and  cloud  water  vapor,  respectively,  measured  in  gm  m"^.  They  are  related  to 
Qc  »  Oh',  q,  by 


Qc  =  f'2,'  ^ 


f  y 

<5s'n^  j 

The  terms  involving  ki  and  kz  are  empirically  derived  for  the  calculation  of 
conversion  and  collection,  respectively.  V  is  the  terminal  velocity  of  liquid 
water  droplets.  It  is  related  to  by 

I  \o.»vr 

\i  _  (20) 


where  K3  is  a  constant. 


,  T  >  k 

M-  ,  T  ^  K 


is  the  density  at  the  base  of  the  cloud. 

We  may  relate  the  horizontal  transport  terms  to  the  entrainment  in  the 
following  heuristic  fashion. 

In  cylindrical  coordinates 

t  /■()),  +c);) 

+  -  (  «v  +  CPs');] 

“STt. 

where  U  is  the  inward  directed  velocity,  and  where  the  subscripts  e  and  i 
refer  to  the  outer  edge  and  the  center  of  the  cloud,  respectively.  Thus 
Ar  =  r  =  radius  of  cloud. 

Since  Qj.g  =  0, 

-  ^  (Qc;  +  [  C  Q,,  -  +  G,,.] 


If  we  put,  for  the  time  being, 
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■,  ~  Mt  (24) 

(These  assumptions  will  be  discussed  later  in  connection  with  entrainment),  then 

-  (w^  +-^iy)(  Qc,'  ~  W-C.-  [cj,.  4  (  0,^.-  a.',]  (25) 

Similarly, 


(  "^y  ^ -5. 1  ^4;  y 


since  Qhe  =  0. 

Then  the  equations  for  and  become 

and  . 

#  TT ‘ j  (28, 

In  this  formula,  we  assume  'X  ^ 

To  derive  the  equation  for  the  vertical  velocity,  we  consider  the  third 
equation  of  motion 

ilt  P  6  6 

Assuming  that  the  environment  is  in  hydrostatic  equilibrium,  that  the  pressure 
gradient  of  the  perturbation  (  can  be  neglected  compared  with  that  of 

the  environment,  and  that  the  drag  forces  include  those  due  to  Qr  and  Qj_  , 
fqudtion  (29)  may  be  written. 

cit  “  (.1+  )T^  -J 

In  Equation  (24),  we  have  (tentatively)  assumed 

Then  we  need  an  equation  for  r.  We  derive  a  continuity  equation  starting 


=  'Tt/l,  P,  lATc. 


(31) 


F 


m 


Then 


fik 

A  -i: 


iT^■'  p.  iAT,.') 


(32) 


[See  Warner  (1970)].  Let  the  average  radial  velocity  around  a  given  cross- 
section  of  cloud  be  U^.  Then  the  average  inflow  is  Upds,  The  actual  inflow  is 

vir 

\  'J-TT/i  Uft  .  (33) 

^  O 

Now  assume  that  the  average  radial  velocity  is  proportional  to  the  vertical 
velocity  in  the  center  of  the  plume. 


laTc-  (34) 

Since  Uf^  -  q.  iaTcj  ,  the  entrainment  is 

"vlT  )T.i  Uft  -  '2-Tr/i;  cv  . 

The  flow  of  mass  into  the  cloud  is  then 

c  vv>"'  sec'* 

Thus  the  convergence  of  mass  flux  is 

=^-IT  fl.q-  p,  (35) 

Using  the  definition  of  in  Equation  (32), carrying  out  the  indicated 

differentiation,  assuming  a/ ,  and  dividing  by  we  find 


If  we  calculate 


__L 

He,-  ^  ^ 


P; 


^  C^,;)  , 


(36) 


we  find 


(37) 


_L  a  Kf 

he 


)w  assume  that 


_ L 

^  ~  Wnc- 


(  dl  He:  _ 

H  c  I  ^ 


Perkey  and  Kreitzberg  (1972)]. 

Equations  (10),  (13),  (20),  (27),  (28),  (30),  (36),  (39)  form  the  system 
the  unknowns  Tc,<?s,V)  Q*  ,(34v,ur,  fx-  solving  the  time  dependent 
>  of  these  equations,  we  write  "f  ^^2 

e  we  have  already  related  horizontal  transport  to  entrainment. 

“ical  Equations 

We  list  these  in  the  order  and  in  the  form  in  which  they  will  be  solved 

Q^.  I  ~  "  *1: 4>i..srfs  1 

--^VO  J 


-♦-‘i.yTS'i 


( ) 'f.'dbev- -»  L  ce  r  wdt«?r  —  uc  [t 

e  n  means  time  step  and  i  means  level.  Here  p  is  pressure  and  ^  .621. 

L  '  '  J  (i 

-  IMT*)  A  1 


'  ti;  4-  q-  4  2  ur  )  4(%p/, ,  -  ^  PcO^ 


?,  =  Jfi 


(53) 


IS  ’ 

-  tu 


A-E 

( o-(>\ 


n:\,\^..^(nir  O 

^1 


-fCV<)<:  -  Ur-  )l,  At 


“1 


A.-Z. 


(54) 


Procedure  for  Solving  S^sjtem 

Given:  or, ,<^r, <^,Cp ,  j  7r , e  ,  •^. ,a  ,-8i.  ,  ,Kj  are  constants 

|<^  r  \S'.^‘1  before  ice  nucleation 
_  U  after  nucleation 

L  ~  'S'«i  r  cal  ^VH-'  ^  T  >  ^  ?r  !<  L,  =1-4  *o 

I  _  '  t>  Ca'  '  7 

•(l,  ,  •80  «!  }k-  V 

initial  Data:  T..  .T,"  .j®,- Y  ,  (  Q^' •)■  , 

Boundary  DaU  Given:  ,  (>;  ,  ^r;  ,  n.'  ,  ( (^X  .  (HiT,  . 

1)  Calculate  from  Equation  (40)  or  (41).  ^ 

2)  Calculate  (  )  fran  Equation  (43).  Use'^=‘*^„at  i-2*  . 

3)  Solve  for 

4)  Calculate  )\T^,  from  Equation  44).  Use^s;^6  at  "Z.  ='a  »,  ,  a'•^a 

are  known  from  Step  3.  f4,  dna»  p*  are  calculated  fron  Equation  (45)  and  (46). 
Stop  calculation  at  level  where  o  . 

5)  Solve  for  /i(  from  Equation  (47). 

-pMl  V  •v/'  —  " 

5)  Calculate  *  ~  '  f  +  Mr^;  ~~_T.  I  from  Equation  (48). 
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Air 


A? 


,  7)  Solve  Equation  (49)  for  .If  Include 

LjtCQti  +  G'«.,y  and  (A  <lsiV  water — »  ice.  Let  L  s  Lj  at  that  level 
and  all  levels  above.  ° 


Vy  ^ 

8)  If  Ti  return  to  Step  6  and  recalculate 

-r  T  *'  “  /  -r  A 

Xi - II  lATc  (Tuj _ -T;  ) 


including  freezing  terms. 


9)  Solve  Equation  (50)  for  (  0,(- ) 


Include  the  term 


If  Y  ^  term. 

In  either  case,  include 

10)  Solve  for  from  Equation  (51). 

If  T  \<  ,  \<^  =  . 

>7-3'  K,  l<i  =  11-^. 

11)  Solve  Equation  (52)  for  (Oa;')  . 

12)  Solve  Equation  (53)  for 

13)  Solve  Equation  (54)  , 

14)  Return  to  Step  1,  but  r\(M  and  IQ^- )  may  not  be  zero. 

15)  Skip  Steps  2,  3,  4,  and  go  to  Step  5. 


16)  Continue  the  procedure  through  Step  13. 

17)  Go  to  Step  14,  and  update  everything  by  1 

culations  at  timeK+\at  the  level  where  _  q  , 


Stop  all  cal- 


Experiments 

We  have  run  eight  experiments,  with  five  sets  of  initial  data.  With 
each  set  of  initial  data,  two  runs  were  made  -  one  with  variable  at,  one 
with  constant  at.  The  variable  at  was  determined  from 

~  (55) 

at  each  time  step.  In  each  case,  probably  due  to  choice  of  initial  conditions, 
the  cloud  life  was  short.  In  an  effort  to  determine  whether  this  brevity  was 
related  to  the  numerical  aspects,  we  ran  each  case  with  several  time  steps, 
each  one  constant  throughout  the  run.  We  selected  the  time  step  which  gave 
the  longest  life  cycle  for  comparison  of  results  using  variable  At.  We  used 
az  =  200  m. 


In  all  but  two  cases,  we  have  chosen  as  boundary  data  w  =  1  m  sec"\ 
r  =  1  km,  Qr  ■  9h  ~  ®  cloud.  We  used  or-  0,15  in  the  entrainment 

equation.  The  initial  data  have  been  varied  in  each  experiment.  The  top  of 
the  cloud  is  that  level  at  which  w  -  0, 

Case  1:  In  an  earlier  paper  (Berkofsky,  1974),  we  have  solved  the  steady- 
state  equations.  We  used  the  results  of  that  experiment  as  input  data  for  T^, 
Tpi  Qrt  Qh*  Fig.  1  shows  the  temperature  profiles  within  the  "cloud"  at  t  «  0^ 
21,  l40,  264  seconds,  for  variable  At.  The  ordinate  is  temperature,  T,  the 
abscissa  the  number  of  the  level,  L.  By  time  264  seconds,  there  are  only  4 
levels  left.  At  the  same  time  (Fig.  2),  the  cloud  radius  at  the  top  has  de¬ 
creased  from  about  1.84  km  at  the  top  to  ,79  km  at  the  top.  The  vertical  vel¬ 
ocity  (Fig.  3)  at  the  last  non-zero  level  has  decreased  from  5.55  m  sec"^  to 
1.87  m  sec-i.  The  cloud  water  (Fig.  4)  has  actually  increased  at  the  top  level 
at  264  seconds,  but  (Fig,  5)  the  hj^rometeor  water  has  decreased  substantially. 
Notice  that  the  patterns  of  r,  w  and  Qc  appear  to  oscillate  widly  at  t  =  140 
seconds,  but  then  settle  down  as  the  cloud  dissipates. 

We  ran  this  case  with  w  at  base  s  105  cm  sec.,  with  similar  results. 

Case  2:  This  is  a  repeat  of  Case  1,  but  with  At  =  9  seconds  -  constant. 

The  cloud  lasts  longer  (until  7  minutes).  The  results.  Fig.  6-10,  for  9,  90, 
and  315  seconds  are  plotted.  These  are  similar  to  Case  1,  but  do  not  oscillate 
as  wildly.  The  increase  in  stability  within  the  cloud,  plus  the  decrease  of 
vertical  velocity  at  the  top,  lead  to  its  destruction,  despite  the  increase  in 
cloud  water. 

Case  3:  We  have  used  real  data  from  Sede  Boqer,  2  March  1982,  at  1322. 

The  environmental  temperature  profile,  Tg,  is  the  same  as  the  curve  labelled 
Initial  in  Fig.  11,  except  that  the  lowest  five  levels  have  been  given  a  temp¬ 
erature  increase  in  order  to  provide  a  perturbation  for  development.  This 
perturbation  was  calculated  from 

:z  \  -i-  cos  ^56) 

where  K  =  1,  2,  3,  4,  5,  (See  Hill,  1974).  ^ 

With  variable  At,  the  cloud  "disappeared"  after  3  minutes.  The  results 
for  18  and  145  seconds  are  shown.  At  that  time,  the  "cloud"  was  still  4400  m 
trdck.  At  the  time  the  data  were  taken  (1322  LST),  there  was  a  cloud  cover  of 
6/8  cumulus  congestus.  But  the  lapse  rate  in  the  lower  layers  become  very 
•■.table,  leading  to  the  clotid's  disappearance.  Neva'theless,  the  cloud  radius, 
the  vertical  velocity,  the  liquid  water,  and  hydrometoor  water  all  reached 
maxiiinum  values  within  the  cloud  at  about  2800  m  before  dissipating  (see  Figs. 
11-15  inclusive).  In  actual  fact,  the  cloudiness  persisted  until  evening. 

This  is  not  predicted  by  the  model. 

Case  4:  Same  as  Case  3,  but  with  At  =  7  seconds  =  constant.  In  this  case, 
the  cloud  lasted  for  3.5  minutes.  The  results  are  shown  for  7  seconds,  20 
seconds,  140  seconds,  245  seconds,  (see  Fig.  16-20  inclusive).  The  results  are 
very  similar  to  those  for  variable  At  (Case  3),  except  they  are  slinntly  more 
stable. 

Case  5;  We  have  used  real  data  for  10  March  1982  for  Sede  Boqer,  at  1257. 
The  calculations  showed  the  cloud  dissipating  after  about  4.5  minutes.  The  re¬ 
sults  for  20,  126,  and  189  seconds  are  plotted  on  Figs.  21-25  inclusive.  Even 
though  the  cloud  had  expanded  at  the  top,  and  the  vertical  velocity  and  liquid 
water  increased,  the  lapse  rate  was  so  stable  that  the  cloud  could  not  grow.  In 
actual  fact,  cloudiness  of  7/8  stratocunulus  persisted  throughout  the  day,  with 
occasional  light  rain.  One  of  the  problems  here  is  that  the  radiosonae  wont 
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through  the  clouds,  giving  cloud  temperature  Instead  of  environmental  temperature. 

Case  6:  Same  ls  Case  5,  but  with  at  *  9  seconds  *  constant.  Figs.  26-30 
Inclusive  summarize  the  results.  The  main  difference  between  this  run  and  that 
for  variable  at  is  that  the  cloud  lasted  about  2  minutes  longer,  and  the  hydro¬ 
meteor  water  stayed  higher  for  a  longer  time. 

Case  7:  We  have  used  real  data  from  Sede  Boqer,  17  March  1982,  at  1300. 

There  was  a  cloud  cover  of  6/8  stratocumulus,  which  disappeared  later  In  the 
day.  With  variable  at,  the  cloud  disappeared  after  3.5  minutes.  The  plots 
Figs.  31-35  inclusive,  show  results  at  t  *  20,  124  seconds.  Actually,  at  124 
seconds,  the  clouds  top  was  still  at  11,200  m,  but  dropped  to  zero  in  another 
1.5  minutes.  In  the  actual  atmosphere,  the  clouds  dissipated  during  the 
afternoon,  but  it  is  not  known  exactly  when.  Notice  that  in  the  calculations, 
the  radius  of  the  top  had  increased  considerably  by  t  «  124  seconds,  the  over¬ 
all  vertical  velocity  had  increased,  as  did  the  top  cloud  water  content.  But 
the  hydrometeor  water  increased  only  in  the  first  few  hundred  meters  above 
the  base,  where  the  cloud  radius  was  smallest  and  where  the  cloud  lapse  rate 
had  become  very  stable. 

Case  8:  Same  as  Case  7,  but  with  at  «  9  seconds  *  constant.  The  cloud 
lasts  for  414  seconds,  instead  of  205  seconds.  However,  the  collapse  of  the 
cloud  took  place  at  an  early  stage,  the  cloud  top  dropping  to  4,000  m  by 
135  sec.  This  is  probably  more  realistic  than  with  variable  at,  since  it  is 
unlikely  that  the  cloud  tops  grew  to  11,200  m  in  the  atmosphere.  The  results 
are  shown  in  Figs.  36-40,  inclusive. 

Inclusion  of  Clwd  Model  In  A  Model  of  Larger-Scale  Flow 

In  an  earlier  paper  (Berkofsl^,  we  have  described  how  a  cloud  model 

might  be  included  in  one  for  the  larger-scale.  In  essence,  we  solve  three 
systems  of  equations:  one  for  the  sub-cloud  layer,  one  for  the  free  air,  one 
for  the  cloud.  In  the  operational  scheme,  we  actually  assume  the  existence  of 
a  cloud  within  each  grid  rectangle.  The  modality  of  the  ensemble  is  to  be  pre¬ 
scribed  from  observations.  Thus  in  the  Central  Pacific,  the  cloud  distribution 
is  bifliodal.  Whenever  the  equations  for  the  larger-scale,  free  air,  predict 
at  least  a  conditionally  unstable  lapse  rate,  it  is  assumed  that  a  cloud  en¬ 
semble  exists  at  that  grid  point.  The  cloud  equations  are  then  solved,  using 
input  data  from  the  free  air  model,  and  assumed  radii  of  bases.  The  conden¬ 
sation  heating  term  in  the  first  law  of  thermodynamics  is  calculated  from 


-L  A. 

At 


(57) 


where  the  bar  is  a  horizontal  average  over  a  grid  spacing.  At  is  the  time  step 
for  the  larger-scale  motions,  f|  Is  the  fraction  of  sky  covered  by  the  ith 
population  of  convective  clouds.  Here  is  the  cloud  liquid  water,  which  is 
calculated  by  the  cloud  model.  The  quantity  Is  calculated  from 

_  At 


At 


At 

In  which  we  assume 
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=  o 


(59) 


A  Qc,- 

At 

This  quantity  represents  the  rate  of  change  of  condensed  water  by  the  clouds. 

This  together  with  the  cloud  heating,  is  fed  into  the  larger-scale  equations, 
thus  influencing  the  entire  surrounding  circulation,  and  hence  the  moisture 
supply  for  subsequent  cloud  calculations. 

Complete  details  of  the  model  and  the  procedure  for  solving  the  system 
are  given  in  Berkofsky  (1974).  The  cloud  model  in  that  paper  differs  in  some 
minor  details  from  the  present  one. 

Conclusions 

We  have  developed  a  one-dimensional  time  depenu^at  "plume"  model  to  pre¬ 
dict  the  growth  of  a  cumulus  cloud.  In  all  of  the  cases  we  studied,  using 
both  real  and  artificial  data,  the  "cloud"  had  a  short  life-time-of  the  order 
of  several  minutes.  This  may  be  due  to  our  choice  of  initial  conditions.  It 
is  possible  that,  had  we  varied  the  initial  temperature  perturbation,  the  initial 
vertical  velocity  at  cloud  base,  the  initial  cloud  radius,  development  would  have 
been  otherwise.  Our  choice  of  a  cloud  radius  and  a  vertical  velocity  at  base 
implies  that  a  cloud  exists.  What  is  more  likely  is  that  the  temperature  diff¬ 
erence  between  our  developing  "cloud"  and  the  real  atmosphere  was  such  as  to  per- 
clude  any  further  development. 

During  the  entire  course  of  the  numerical  experimentation,  we  never 
reached  the  point  where  we  had  sufficient  observations,  both  in  the  ambient 
atmosphere  and  of  the  actual  cloud  growth,  to  test  the  predictions.  Even  on 
days  when  there  were  clouds,  there  were  no  continuous  cloud  observations,  so 
it  is  not  certain  whether  a  specific  cloud,  predicted  to  decay,  actually  did  so. 

Nevertheless,  it  is  feH  that  the  model  deserves  further  testing.  In  the 
very  near  future,  we  will  be  in  a  much  better  position  to  measure  required  data 
as  input  and  verification  of  the  model.  These  will  include  radar  data. 
Recommendations 

It  is  recommended  that  tiie  model  be  re- tested  when  more  data  become  avail¬ 
able  in  the  Negev.  It  is  possible  that  such  data  exist  in  U.S.  deserts. 

It  is  also  recommended  that  the  model  be  used  in  conjunction  with  a 
larger-scale  model.  Thus,  when  the  larger-scale  model  predicts  that  the  air 
at  a  grid-point  is  at  or  near  saturation,  that  the  lapse  rate  is  conditionally 
unstable,  that  the  vertical  velocity  of  the  larger-scale  motion  is  upward,  it 
is  assumed  that  an  ensemble  of  clouds  exist.  The  modality  and  the  input  data 
(radius  of  base,  vertical  velocity  of  base,  cloud  water  and  liquid  water  at 
base)  will  be  assumed  from  known  cloud  statistics  of  the  region,  and  will  be 
used  as  input  data,  together  with  ambient  temperature  from  the  larger-scale 
model,  to  solve  the  cloud  equations.  The  total  latent  heat  released  at  each 
grid  point  is  then  calculate  from  Eq.  (58),  and  fed  into  the  first  law  of 
thermodynamics  for  the  larger-scale  motions. 

Finally,  it  is  recommended  that  the  model  be  modified  to  include  a  dust 
equation.  As  a  first  attempt,  we  might  consider  a  dry  atmosphere,  lacking 
either  cloud  water  or  hydrometeor  water,  in  which  dust  is  entrained  during 
conditionally  unstable  conditions.  It  would  be  of  great  interest  to  calcul¬ 
ate  the  height  to  which  a  dust  cloud,  emanating  from  the  ground,  would  grow. 

Even  a  steady-state  calculation  would  be  revealing. 
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Content  x  10  (gmgm-1)  vs.  Level  Number  (Vertical  Resolution 
a  Function  of  Time  (seconds).  Case  1. 


200  m)  as  a  Function  of  Time  (seconds).  Case  2. 


eioc 


Cloud  Water  Content  xlO  (gnqm"^)  vs.  Level  Number  (Vertical  Resolution 
is  200  m)  as  a  Function  of  Time  (seconds).  Case  2. 


T«18 


Fig.  11.  Cloud  Temperature  xl0“2  (Degrees  K)  vs.  Level  Number  (Vertical  Resoluti 
is  200  m)  as  a  Function  of  Time  (seconds).  Case  3. 


Fig.  14.  Cloud  Water  Content  xlO^  (gnigni-l)  vs.  Level  Number  (Vertical 
Resolution  is  200  m)  as  a  Function  of  Time  (seconds).  Case  3. 


Resolution  is  200  m)  as  a  Function  of  Time  (secon 


Temperature  xlO"^  (Degrees  K)  vs.  Level  Number  (Vertical  Resolution 
I  tn)  as  a  Function  of  Time  (seconds).  Case  4. 


•city  xlO-^  (cm  sec"^)  vs  Level  Number  (Vertical  Resolution 
a  Function  of  Time  (seconds).  Case  5. 


Fig.  28.  Vertical  Velocity  xlO-2  (cm  sec"M  vs.  Level  Number  (Vertical  Resolution 
Is  200  m)  as  a  Function  of  Time  (seconds).  Case  6. 


Hydronieteor  Water  Content  xlO'^  (gmgm"^)  vs.  Level  Number  (Vertical 
Resolution  is  200  m)  as  a  Function  of  Time  (seconds).  Case  6. 


I  Radius  xlO"^  (cm)  vs.  Level  Number  (Vertical  Resolution 
Function  of  Time  (seconds).  Case  7. 


:esolution  is  200  m)  as  a  Function  of  Time  (seconds).  Case 


Temperature  xlO"^  (Degrees  K)  vs.  Level  Number  (Vertical  Resolution 
I  m)  as  a  Function  of  Time  (seconds).  Case  8, 


'elocity  xlO"^  (cm  sec~l)  vs.  Level  Number  (Vertical  Resolution 
as  a  Function  of  Time  (seconds).  Case  8- 


Appendix  -  Computer  Program 

The  attached  program  is  in  two  parts.  PROGRAM  CLOUD  1  goes  through  the 
initialization  procedure  to  calculate  the  initial  values,  and  copies  the  re¬ 
sults  to  files.  PROGRAM  CLOUD  2  carries  out  the  calculations,  according  to 
the  steps  outlined  in  Procedure  for  Solving  the  System. 

The  calculations  were  done  on  the  PDP  1134  computer  at  Kibbutz  Sede 
Boqer.  The  program  was  written  in  BASIC  PLUS  language. 
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EXTEND 

I - !  j 

i  PROGRAM  CLOUDl  IS 

! - 1 

DIM  TE(61X>rQE(61X)fT<61Zr2X)»QC<61Xf2X)fQH(61%r2X>  t 
DIM  aS(61X)rU(^lXf2X)rR(61Xr2X>fR0(61Xr2X)fP(61Xf2X>  t 
DIM  MU(61Xr2X>»V(61Xr2X)»aSI(61X)fDQS(61X) 

I - !  I 

»  CONSTANTS  AND  DEFINITIONS  It 

j - 1 

ALFA-0.15  \  Al-2.39*<10*»<-8))  \  G-980  \  CP-0.239  t 

J-4.185*<10*#7)  \  RC-2. 87* <10*#6)  \  EPS-0.421  t 

A»0.1*(10**(-5)) 

, - 1  j 

I  INPUT  DATA  FROM  TERMINAL  I  t 

j - j 

PRINT  'CLOUDl  •>TIME*<OX)>'  'JDATE^COX)  \  PRINT  S 


INPUT  'DZ  VALUE' »DZ  \  PRINT  t 

INPUT  'LEVELS  NUMBER  -  MAXIMUM  60'»LEVX  \  PRINT  t 
INPUT  'U<1»1)  VALUE'»W(1X»1X)  \  PRINT  t 
INPUT  'P<1»1)  VALUE' ?P<1X» IX)  \  PRINT  S 
INPUT  'FILE  NUMBER' »D*  \  PRINT 


, - !  i 

!  PRINT  DATA  IS 

I - ! 


PRINT  \  PRINT  '  CONSTANT  S'  S 

■  ssssmtssBBsaBas— — BsasBsssssrassss  *  \  PRINT  S 

PRINT  '  ALFA»')ALFA)'  A1»')A1)'  G-'>G>'  CP-')CP»'  J='»J»  S 

'  RC-')RC?'  EPS-'?EPS»'  A-')A  \  PRINT  S 

PRINT  '  DATA  FROM  TERMINAL'S 

”  ssaBaaevsnaBaaBarsESBaBsasssasstsflssBBasssassasBSSSsssBBsscssssESSsasssscsss:  ”  \  PRINT  & 

PRINT  '  W(l»l>»')W(lXrlX)) '  P( 1 f 1 ) - ' )P< IXr IX) » '  D2-'»D2)  S 

'  LEVELS-' »LEVX»'  FILE  NUMBER-' )D*  \  PRINT  \  PRINT 


! - 1  4 

I  OPEN  FILES  AND  COPY  DATA  TO  ARRAYS  !  S 

, - ( 


OPEN  'CLTE.'+D*  FOR  INPUT  AS  FILE  IXfMODE  8192X  S 
INPUT  *lXfTE(IX)  FOR  IX-IX  TO  LEVX  \  CLOSE  IX  S 
OPEN  'CLQE.'+D*  FOR  INPUT  AS  FILE  IXfMODE  ei92X  S 
INPUT  *lXfQE<IX)  FOR  IX-IX  TO  LEVX  \  CLOSE  IX  S 
OPEN  'CLT.'+D*  FOR  INPUT  AS  FILE  IXfMODE  ei92X  S 
INPUT  *lXfT(IXflX)  FOR  IX-IX  TO  LEVX  \  CLOSE  IX 

, - !  I 

I  PRINT  DATA  FROM  FILES  IS 

I - 1 

STl*-'  LEVEL  '  \  ST1*-ST1*+'  TE<I)  QE<I)  '  S 

8T1*»8T1*+'  T(IfN)'  \  ST2*-'  -  '  S 

8T2*-ST2*+' -  '  FOR  IX-IX  TO  3X  S 

ST3*-'  *♦'  \  ST3*-ST3*+'  FOR  IX-IX  TO  3X  \  PRINT  S 

PRINT  '  DATA  FROM  FILE  S'  S 

PRINT  *  BaBaiSBaiHinB»aBttfliatBB8l88BSaiflBM9SSBI«a8«aBflBBn8B8BSB«aiB  *  \  PRINT  1 

PRINT  STl*  \  PRINT  8T2*  S 

PRINT  USING  ST3*fIXfTE<IX)fQE<IX)fT<IXflX)  FOR  IX-IX  TO  LEVX 


, - !  4 

I  CALCULATE  INITIAL  VALUES  I  S 

I - 1 


T(lXf2X)-T(lXflX)  \  QC(lXflX)fQC(lXf2X)-0  S 
QH(lXflX)faH(lXf2X)-0  S 


505 

\ 


\ 

550 


555 

\ 

\ 

\ 

\ 

<&00 


<S05 

\ 

\ 

650 


655 

\ 

\ 

700 


710 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 


790 

\ 

\ 

820 

1000 


1005 

\ 

\ 

\ 

\ 

N 

\ 

\ 

\ 

\ 

\ 


RdZf  IZ)  rR(lZ»2Z)  =  10*!|t5  \  P<  IZr  2Z)  »P  <  IZ » IZ) 

, - f  i 

I  STEP  1  -  CALCULATE  P(I»N)  FOR  N»l»  FOR  1-2  TO  LEO  !  S 

I - ! 


FOR  IZ-2Z  TO  LEVZ  t 
EZl-0  I 

EZ1»EZ1+DZ/T<KZ»1Z)  FOR  KZ-IZ  TO  IZ  I 
P(IZ»1Z)»P<1Z»1Z)*EXP<  <-G/RC)*EZl)  t 
NEXT  IZ 


! - !  s 

!  STEP  2  -  CALCULATE  RO  FOR  N-l»  FOR  1=1  TO  LEV  !  S 

I - ! 

RO<IZ»lZ)«P<IZ»lZ)/<RC»T(IZrlZ) )  FOR  IZ=1Z  TO  LEVZ  S 
RO ( LEVZ-t- 1 Z » 1 Z  )  -2«R0 ( LE VZ  r  IZ )  -RO  ( LEVZ- IZ » IZ )  S 
R0<1Z»2Z)-P<1Z»2Z)/<RC*T(1Z»2Z)) 

I - !  i 

I  STEP  2A  -  CALCULATE  QCfQHrV  FOR  N-l»  FOR  1=1  TO  LEV+1  !  t 

j - ! 

QC(IZflZ)»0  FOR  IZ-IZ  TO  LEVZ+IZ  t 
aH(IZ»lZ>=0  FOR  IZ-IZ  TO  LEVZ+IZ  S 
V(IZrlZ)-0  FOR  IZ-IZ  TO  LEVZ+IZ 

I - !  I 

!  STEP  3  -  CALCULATE  MU  -  EQUATION  A7  !  t 

!  STEP  5  -  CALCULATE  U  -  EQUATION  43  !  S 

!  STEP  5A  -  CALCULATE  R  -  EQUATION  44  !  S 

!  STEP  4  -  CALCULATE  QS  -  EQUATION  40  !  i 

!  STOP  CALCULATION  WHERE  U-<0  I  t 

I - ! 


FOR  IZ-IZ  TO  LEVZ  t 
MZ-IZ+IZ  i 

MU(IZ»lZ)=2*ALFA/R(IZf IZ)  \  E2-L0G10(T< IZ» IZ) )  t 
El— 2937 . 4/T <  IZ » IZ )  -4 . 9283*E2+22 .5518  & 

E3-P(IZdZ)*0.0001  t 
QS(IZ)-(EPS/E3)«(10««E1)  t 

El-<l+0.61*QS<IZ))*T<IZf IZ)  \  E2-<1+0.61*QE(IZ) )*TE(IZ)  i 
E3-2*MU<IZf  1Z)#<W<IZ»1Z)**2)JKDZ  \  E4=(E1-E2)/E2  S 
E5-(U<IZrlZ)««2)F2«G«0Z«E4-E3  i 
IF  E5>-0  AND  W(IZrlZ)>0  t 
THEN  790  t 

ELSE  LEVZ-IZ  \  GO  TO  1000 

W(HZf 1Z)-SQR<E5)  \  E1-L00(U(HX» IZ ) )-LOG(U( IZf IZ) )  S 
E2»L0G<R0<MZ»1Z))-L0G(R0<IZ»1Z))  \  E3=R( IZ» lZ)/2  S 
R<MZ»lZ)-R<IZrlZ) +ALFA*DZ-E3  <  E 1 +E2 ) 


NEXT  IZ 

I - 

!  PRINT  RESULTS 

I - 

STl*-*  LEVEL  TE(I)  ■  * 

ST1*-ST1*+*  QE(I)  T(JfN)  •  t 

8T1$-8T1*+*  QCdfN)  QH<I»N)  *  t 

8T1*-8T1*+*  QSdrN)  WdfN)  •  \  ST2t-* - 

8T2t-8T2f+* -  •  FOR  IZ-IZ  TO  7Z  \  ST3*-*  t 

ST3*-ST3f+*  FOR  IZ-IZ  TO  7Z  t 

ST4*«*  LEVEL  RdrN)  '  » 

ST4f-ST4*+*  ROdfN)  PCIrN)  *  t 

8T4f-8T4*+*  MUdfN)  VdrN)  •  t 

8T4f-8T4»-l-*  Q81d»N)  DQ8drN)  *  \  PRINT  t 

PRINT  *  RE8ULT8  -  INITIAL  VALUE8*  t 


I  t 
!  t 

I 


t 


nn Ayoo r 


1 


\  PRINT  *  aaasBaiaBBBKBaaasaaBKaBaaaBassBaesasssaBBB  *  K 

\  PRINT  \  PRINT  STl*  \  PRINT  ST2*  t 

\  PRINT  USING  8T3«rIX»TE<IX>rQE(lX)rT(IX»lX)rQC<IXflX)fQH(IXflX)»  S 

QS(IX>fW(IX»lX>  FOP  ix-ix  TO  LEUX  S 
\  PRINT  \  PRINT  \  prim  \  PRINT  ST4*  \  PRINT  ST2*  t 

\  PRINT  USING  ST3«rIXrR(IXrlX>rR0(IX»lX>FP(IXflX)fHU(IXf lX>f  I 

U(IXflX)»QSI(IX)fOQS(IX)  FOR  IX«1X  TO  LEUX  S 
\  PRINT  \  PRINT 


1100  I - !  t 

I  COPY  RESULTS  TO  FILES  I  t 

f - j 


1105  INPUT  'SAVE  RESULTS  IN  EXTERNAL  FILES  <YES  OR  NO>’>DD»  t 

\  GO  TO  10000  IF  DD*»‘NO'  t 

\  GO  TO  1120  IF  DD*-'YES*  t 

\  PRINT  ‘INCORECT  ANSWER*  \  GO  TO  1105 

1120  OPEN  'CLP.'+D*  FOR  OUTPUT  AS  FILE  IX  t 

\  PRINT  OlXfPdXflX)  FOR  IX-IX  TO  LEVX  \  CLOSE  IX  t 

\  OPEN  'CLOS^'+D*  FOR  OUTPUT  AS  FILE  IX  i 

\  PRINT  OlXrOSdX)  FOR  IX-IX  TO  LEVX  \  CLOSE  IX  S 

\  OPEN  'CLW.'+D*  FOR  OUTPUT  AS  FILE  IX  t 

\  PRINT  #lXfWdXflX)  FOR  IX-IX  TO  LEVX  \  CLOSE  IX  > 

\  OPEN  'CLR.'+D*  FOR  OUTPUT  AS  FILE  IX  t 

\  PRINT  #lXrRdX»lX)  FOR  IX-IX  TO  LEVX  \  CLOSE  IX  t 

\  OPEN  ■CLMU.*+D*  FOR  OUTPUT  AS  FILE  IX  t 

\  PRINT  ♦IXfHUdXrlX)  FOR  IX-IX  TO  LEVX  \  CLOSE  IX  i 

\  OPEN  'CLRO^'+D*  FOR  OUTPUT  AS  FILE  IX  i 

\  PRINT  flXfROdXdX)  FOR  IX-IX  TO  LEVX  \  CLOSE  IX 

10000  END 


wiMiidMfiiiiiiifi 


■liMlMMIliill 


1  EXTEND 

5  ! - !  X 

!  PROGRAM  CL0UD2  !8 

I - 1 

10  t - !  & 

!  DEFINITIONS  !S 


15  DIM  TE(61X)f  QE(61Z)»  T<61%»2Z)»  QC<61Z,2X)f  QH(61Zf2%)  S 

\  DIM  QS(61Z)»  W(61Zf2X)»  R(61%f2%)»  R0<61X»2Z)»  P(61Z»2%)  S 

\  DIM  MU(61%»2Z) »  V(6iXf2X)f  QSI(61X)»  DQS<61%) 

50  ! - !  S 

!  CONSTANTS  !& 


55  ALFA=0.15  \  Al=2.39»( 10**(-8) )  \  G=980  \  CP=0.239  S 

\  J=4.185»<10»*7)  \  RC=2.87*<10**6)  \  EPS=0.621  & 

\  A=0. 000001  \  E78=7/8  % 

\  ST1*='  LEVEL  T(I»N)  W(I»N)  R(IfN)  ’  S 

\  ST1*=ST1*+’  P(I,N)  QC(I»N)  GH(I»N>  '  S 

\  ST2*=*  LEVEL  RO(I»N)  MU(IfN)  V(I»N)  '  & 

\  ST2*=ST2$f  QS<I)  QSI<I)  DQSd)  ■  « 

\  ST3*=‘  -  -  -  -  ■  & 

\  ST3$=^ST3$’f "  _ _ "  & 

\  ST4$=*  ♦♦  ■  S 

\  ST4*=ST4$f  ■  Z 

\  SUMTIM=0 

100  ! - !  Z 

!  INPUT  DATA  FROM  TERMINAL!* 

1 - I 


105  PRINT  ’C  L  0  U  D  2  '  rTIME$(OZ)  J  ■■  •fDATE$(OZ)  \  PRINT  & 

\  INPUT  ’DZ  VALUE" ?DZ  \  PRINT  * 

\  INPUT  "LEVELS  NUMBER  -  MAXIMUM  60"JLEVZ  \  PRINT  * 

\  INPUT  "TIME  STEPS  -  MAXIMUM  250"fSTEPS%  \  PRINT  X 

\  INPUT  "FILE  NUMBER" fD$  \  PRINT  * 

\  INPUT  "CHANGE  T(IfN)  <YES>  >  <NO>*?DD$  \  PRINT  * 

\  IF  DD$="NO"  Z 

THEN  110  X 

ELSE  INPUT  "RANGE  TO  CHANGE" »D 
110  INPUT  "IS  DT  CONSTANT  <YES>  f  <NO>"»DDT$  Z 

\  IF  DDT$="NO"  Z 

THEN  200  Z 

ELSE  INPUT  "DT  VALUE" fDT 


200  ! - !  * 

IPRINTDATA  !X 


205  PRINT  \  PRINT  "  CONSTANTS"* 

\  PRINT  •  ================:====.==:===«=•  \  print  Z 

\  PRINT  "  ALFA=" JALFA" f  Al="fAlf"  G="fGJ"  CP="fCPf"  J="fJfX 

"  RC="fRCf"  EPS="fEPS>"  A="JA>"  TIME= " ? SUMTIM  Z 
\  PRINT  \  PRINT  "  DATA  FROM  TERMINAL" 

\  PRINT  ■  =  =  =  ====  =  =  ==  =  =  =  *=====  =  =  =  ====  =  =  =  =  =  ==  =  =  =  =  =  =  =  ==  =  ==  =  =.•=.=.•==  =  =:=: 

\  PRINT  Z 

\  PRINT  "  DZ="?DZ?"  LEVELS="JLEVZ> •  STEPS=" f STEPSZf  8 


\  PRINT  ■  FILE  NUMBER=">D*  \  PRINT  \  PRINT 

250  ! - !  X 

!  OPENES  FILES  AND  COPY  DATA  FROM  FILES  TO  ARRAYS  !  X 


255  OPEN  "CLTE.'+D*  FOR  INPUT  AS  FILE  IZfMODE  8192Z  X 
N  INPUT  tlZ^TEdZ)  FOR  IZ=1Z  TO  LEVZ  \  CLOSE  IZ  X 


OHAiycc 


OPEN  'CLQE.'+D*  FOR  INPUT  AS  FILE  IZrMODE  8192%  & 
INPUT  #1%»QE(I%)  FOR  I%=1%  TO  LEU%  \  CLOSE  1%  S 
OPEN  •CLP.'+D*  FOR  INPUT  AS  FILE  l%fMODE  8192%  S 
INPUT  #l%»P(I%fl%)  FOR  I%=1%  TO  LEU%  \  CLOSE  1%  i 
OPEN  'CLT.'+D*  FOR  INPUT  AS  FILE  1%»M0DE  8192%  S 
INPUT  ♦1%»T(I%»1%)  FOR  I%=1%  TO  LEV%  \  CLOSE  1%  i 
OPEN  'CLOS.’+Df  FOR  INPUT  AS  FILE  1%»M0DE  8192%  i 
INPUT  ♦l%fQS(I%)  FOR  I%=1%  TO  LEV%  \  CLOSE  1%  S 
OPEN  'CLW.'+D*  FOR  INPUT  AS  FILE  l%fMODE  8192%  S 
INPUT  ♦1%»W<I%»1%)  FOR  I%=1%  TO  LEV%  \  CLOSE  1%  i 
OPEN  'CLR.'+D*  FOR  INPUT  AS  FILE  1%»M0DE  8192%  S 
INPUT  #l%fR<I%fl%)  FOR  I%=1%  TO  LEU%  \  CLOSE  1%  i 
OPEN  'CLMU.'+Df  FOR  INPUT  AS  FILE  l%rMODE  8192%  i 
INPUT  #1%»MU(I%»1%)  FOR  I%=1%  TO  LEV%  \  CLOSE  1%  S 
OPEN  ■CLRO«'+D$  FOR  INPUT  AS  FILE  l%fMODE  8192%  S 
INPUT  #l%»RO(I%fl%)  FOR  I%=1%  TO  LEU%  \  CLOSE  1% 


ZEROS 


AND 


UPDATE 


T(I>N) 


MAT  aC=ZER  \  MAT  QH=ZER  \  MAT  U=ZER  I 
MAT  QSI=ZER  \  MAT  DQS=ZER  i 
IF  DD»='NO*  S 
THEN  510  S 

ELSE  T(I%»1%)=T(I%»1%)+D  FOR  I%=1%  TO  LEU% 


MAIN  LOOP  -  STEPS  OF  TIME 


FOR  T%=1%  TO  STEPS%  !  T%  -  TIME  INDEX  & 
DEG%=0%  \  E1=0 

FOR  I%=1%  TO  LEW.  !  CALCULATE  DT  S 
IF  W(I%»1%)>E1  THEN  E1=W<I%»1%) 

NEXT  1% 

DT=DZ/E1  IF  DDT*='N0‘ 

SUMTIM=SUMTIM+DT 

FOR  I%=1%  TO  LEV%  !  IX  ~  LEVELS  INDEX  i 

IN%=I%+1%  \  IM%=I%-1%  \  IF  T%=1%  THEN  580 


CALCULATE  QS  WATER  EQ.  40 


El=L0G10<T<I%fl%) )  S 

E2=-2937,4/T(I%»l%)-4.9283!|cEl+22,5518  % 
E3«0.0001#P(I%»1%)  S 
QS(I%)=<EPS/E3)*(10##E2) 


CALCULATE  QS  ICE 


E  Q.  41 


El*-2667/T(I%»l%)+9.5553  S 
E3“0.0001*P(I%f 1%)  S 
QSI(I%)  =  (EPS/E3)i|c<10*»El) 


CALCULATE  DQS 


E  Q.  42 


DaS<I%)=QS(I%)-QSI<I%)  S 
IF  I%«1%  S 
THEN  590  % 

ELSE  600 

T<l%f2%)*T(l%»l%)  S 
P<lXf2%)=P(l%»l%>  8 


\ 

\ 

\ 

\ 

\ 

600 


601 


605 

\ 

610 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

N 

650 


655 

\ 

660 


665 
\ 

\ 

666 
667 
\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

669 

700 


R<l%>2%)=R(l%f IZ)  S 
R0(1Z»2%)=P(1%»2Z>/(RC*T(1%>2%)  )  S 

MU<1%»2Z)=2*ALFA/R(1%»2Z)  & 

GO  TO  2000  !  TO  NEXT  IX 


CALCULATE  T  -  EQ.  49 


% 


IF  IZOLEVZ  S 
THEN  605  & 

ELSE  IF  T(IZ»1Z).>T<IMZ»1Z)  S 

THEN  T(INZ»lZ)=2*T(IZf 1Z)-T(IMZ»1Z)  % 

ELSE  T(INZ»1Z)=0*9*T(IZ»1Z) 

LA=677  \  LF=80  \  LS=677  % 

IF  T(IZ»1Z)>273  THEN  LA=595  \  LF=0  \  LS=0 

El=T<INZrlZ)-T<IZ»lZ>  S 

E2=(  (TdNZf  1Z)+T(IMZ>1Z)  )/2) -U<  IZ»  IZ ) * ( DT/DZ )  *E1  S 
E3=W(IZ»1Z)*DT  S 

E4=(Al*G/CP)»(l+QS(IZ)*LA*J)/(RC*T(IZf IZ) )  S 
E5=MU(IZ»1Z)»(T(IZ»1Z)-TE(IZ))  8 
E6=MU(  IZdZ)*  (LA/CP) *(QS ( IZ ) -QE ( IZ ) )  S 

E7=l+( (EPS*(LA**2)#QS(IZ) )/(CP*Al*RC*(T<IZf 1Z)**2) ) )  8 
E8=LF#<GJC(IZ»lZ)+QH(IZf  IZ)  )  8 
E9=LS#DQS(IZ)  8 
E10=1/CP  8 

T  <  I Z » 2Z )  =E2-E3)|c  <  E4+E5+E6 ) /E7"E1 0*  ( E8+E9 ) /E7 


CALCULATE  QC  -  EQ.50 


8 

8 


Kl=0.0015  \  K2=0,0696  8 

IF  T<IZflZ)>273  8 

THEN  K1=0. 00075  \  K2=0.0052 

IF  IZOLEUZ  8 

THEN  665  8 

ELSE  IF  QC(IZ»lZ)>aC<IMZ»lZ)  8 

THEN  QC(INZFlZ)=2*QC(IZf lZ)-QC(IMZf IZ)  8 
ELSE  QC<INZf1Z)=0.9*QC(IZf1Z) 

E1-(QC(INZ»1Z)+0C(IMZf1Z) )/2  8 
FOR  XZ=IMZ  TO  INZ  8 
IF  E1<QC(XZf1Z)  8 
THEN  El=aC<XZFlZ) 

NEXT  XZ 

El*El-U(IZFlZ)*(DT/DZ)*(QC<INZflZ)-QC<IZ»lZ) )  8 
E2»G#W(IZf1Z)*QS<IZ)»DT/(RC#T(IZ»1Z) )  8 
E3=EPS#LA*J#QS(IZ)#DT/(RC#(T(IZf1Z)#*2) )  8 

E4*<T(IZ»2Z)-T(IZ»lZ))/DT+(U(IZFlZ)/DZ)*(T(INZrlZ)-T(IZFlZ) )  8 
E5»MU(IZf l)*W<IZrlZ)*DT»(QS(IZ)-QE<IZ)+QC(IZf IZ) )  8 
E9*R0(lZf lZ)/R0(IZf IZ)  8 

E7®K2#DT*(E9*»0.5)*(R0(IZrlZ)»»E78)*QC(IZFlZ)*(QH(IZF 1Z>**E78)  8 

QC(IZf2Z)=E1-E2-E3*E4-E5-E7  8 

E8=Kl*DT»<QC(IZf 1Z)-(A/R0(IZf1Z) ) )  8 

IF  QCdZf  1Z)><A/R0(IZf1Z)  )  8 

THEN  QCdZF2Z)=QCdZF2Z)-E8 

IF  ClCdZ»2Z><0  8 

THEN  QCdZ»2Z)“0 


CALCULATE  QH  -  EQ.  52 


8 

8 


0n>177C  A 


I%=LEV%  S 

:n  710  I 

IE  725 

:N  QH<IN%>lZ)*2*QH(I%f l%)-aH(IM%rlX)  X 
5E  QH(INZ»1Z)*0.9»QH(I%»1Z) 

R0(IZ»1Z)>R0<IHZ»1%)  S 

:N  R0(INZ»l%)=2*R0(I%f 1Z)-R0<IMZ»1%)  X 
3E  ROdNZf  1Z)*0,9»R0<I%»1Z) 

V<I5:»lZ).>V(IMZrlZ)  S 
iN  V<IN%»lZ)=2*V<IZf l%)-y(IM%f 1%)  X 
3E  V(IN%»l%)=0.9*0<IZrlX) 

3=W(I%»lZ)+V<IXfl%)  X 
=<QH(INXfiX)+QH(IhXflX)>/2  X 
=El-E10*(DT/DZ)*(QH(IN%»l%)-QH(IZd%) )  X 
=V(IZf 1%)*(L0G(R0(IN%»1%) )-LOG<RO(I%fl%) ) )/DZ  X 
=(0(INZfl%)-y(IZ»lZ))/DZ  X 
=MU(IZ»lZ)#W(I%,lZ)*QH(I%f IZ)  X 
=Kl»(QC(IZ»l%)-(A/RO<IZfl%)))  X 
=RO(l%»l%)/RO<IZrlZ)  X 

=K2*(E8**0.5)*(R0(I%»1%)!K*E78>*(QH(I%»1Z>»*E78)  X 
< I%»2%)=El+DT*(-QH(IXf 1%)*(E3+E4)-E5+E9)  X 
QC(I%,1X)>(A/R0<l%f 1%) )  X 
EN  QH(I%»2%)=QH<I%»2Z)+(DT*E6) 

QH(I%»2%)<0  X 
lEN  QH<I%»2%)=0 


CALCULATE  V-EQ.  51 


X 

X 


TdZf  IZ) <=273  X 
lEN  K3=15.39  X 
■SE  K3=11.58 

■  aH(iy.t2%)<=0  X 
lEN  Ud%»2Z)=0  X 

■  SE  Ud%r2Z)=K3*(QHdZf2%)»*0.125)*(-130) 


CALCULATE  W-EQ.  53 


X 

X 


•  IX-LEVX  X 

lEN  IF  Ud%f  l%>>UdMXd%>  X 

THEN  WdN%»lZ)=2*WdZfl%)-WdM%f  1%)  X 
ELSE  WdN%»l%>=WdZ>lZ)/2+l 
,=Wd%f  lZ)*<WdN%»lZ)-Wd%f  1%)  )»DT/DZ  X 
l=d+0.61#QSdZ)  )*TdZ»lZ)  X 
(=(l+0.61*QEdZ)  )*TEd%)  X 
l=QC<IZdZ)+QHdZ»lZ)  X 
i=MUdZf  lZ)*<WdZf  1Z)»*2)*DT  X 
»=<E2-E3)/E3  X 

IZ»2Z)  =  <WdNZf  lZ)+WdMZdZ)  )/2-El+DT*G*(E6-E4)-E5  X 
Ud%f2Z)>0  X 
lEN  820  X 

.SE  LEVZ-IZ-IZ  \  DEGZ»1Z 


CALCULATE  R  -  EO.  54 


X 

X 


I%«LEU%  X 

tEN  IF  Rd%dZ)>RdM%rl%>  X 

THEN  RdN%dX)«2*RdZ»lZ>--RdMZ»lZ)  X 


\0 


ELSE  R(lN%fl%)=0,9*R(r/.fl‘/.) 

El=W<I%f  1X)»DT»(R(IN%»1%)-R(I%»1%)  )/riZ  S 
E2=ALFA*DT*W<IXf IX)  S 
E3=U(IXrlX)*R(IX»lX)#DT/2  S 
E4=(L0C(W(INX»1X) )-LOG(U(IX» IX) ) )/DZ  & 
E5=(L0G(R0(INX»1X) )-LOG(RO(IX»lX) ) )/DZ  & 
R(IX»2X)=(R( INX» 1X)+R( IMX» IX) )/2-El+E2-E3»<E4+E5) 


>0  j - !  x 

!  CALCULATEP-EQ,  46  !S 


>5  E1=-G/RC  \  E2=0  & 

FOR  JX=1X  TO  IX  S 
IF  T(JX»2X)<>0  S 
THEN  E2=E2+DZ/T( JX»2X) 

0  NEXT  JX 

.5  E3=E1*E2  & 

P  ( IX » 2X )  =P  ( 1%  r  2X )  JKEXP  ( E3 ) 


iO  , - I  s 

!  CALCULATE  ROEQ.  45  !& 

i5  R0(IX»2X)=P<IX»2X)/(RC*T(IX»2X) ) 

,0  I - ! 

!  CALCULATE  MU-EQ.  1  !& 

.5  MU(IXf2X)=2»ALFA/R(IX»2X)  S 

IF  DEGX=1X  S 
THEN  2001  S 


too  NEXT  IX  !  END  LEVELS  LOOP  ! 

>01  ! - !  X 

!  PRINT  RESULTS  !X 

I - I 


>02  PRINT  •  OUTPUT  -  STEP  OF  TIME  'fTX?'  TIME= ' J SUMTIM f  S 

■  LEVELS=* fLEVX  & 

IF  TX=1X  OR  TX-(TX/5X)*5X=0X  i 
THEN  2005  i. 

ELSE  2200 
'05  PRINT  S 

PRINT  •  OUTPUT  -•  STEP  OF  TIME  ■fTXf*  TIME^ '  JSUMTIMr  & 

‘  LEVELS='?LEVX  S 

00  FOR  JX==1X  TO  2X  i 

PRINT  S 
IF  JX=1X  S 
THEN  PRINT  ST1$  S 
ELSE  PRINT  ST2* 

15  PRINT  ST3$  % 

IF  JX=2X  i 
THEN  2140 

20  FOR  LX=1X  TO  LEVX  i 

PRINT  USING  ST4*^LX»T(LXf2X) »W(LX»2X) fR<L%r2X) f * 

P(LX»2X)fQC(LXf2X)fQH(LX»2X)  S 

NEXT  LX  8 

GO  TO  2150 

40  FOR  LX*1%  TO  LEVX  8 

PRINT  USING  ST4*fLX»R0(LXr2X)fMU(L%f2X)fV(LXF2X)»  8 


on>i7-7co 


QB(L%)»aSI(LX>»DQS(L%)  S 

\  NEXT  L%  \  PRINT  \  PRINT 

2150  NEXT  J% 


1  -- 

1 

1  — .  - 

TRANSFER 

FROM  COLUMN  2 

TO  COLUMN  1 

f 

IF 

LEVZ=1Z  THEN 

10000  ! 

STOP 

RUN 

2215  FOR  JZ=1X  TO  LEVX  % 

\  r(jXfiy.)=T<JXf2y.)  \  qc(  jz»iz)=qc(  jx»2Z)  \  oh(  jz»iz)=qh(  jz»2Z) 

\  W(  JZ»1Z)=U(JZ»2Z)  \  R<  JZf  1Z)=R(  JZ»2Z)  \  P ( JZ >  IZ )  ==P ( JZ » 2Z )  & 

\  MU( JZ»1Z)=MU< JZ»2Z)\  R0( JZ»1Z)=R0(JZ»2Z)  \  V( JZ » IZ) JZf 2Z )  8 

\  NEXT  JZ 

3000  NEXT  TZ  !  END  TIME  LOOP  ! 

10000  END 
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